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ABSTRACT 

The  Fourier-ray  method  involves  ray  tracing  in  a  Fourier-transform  domain.  The  ray  solutions  are  then 
Fourier  synthesized  to  produce  a  spatial  solution.  Here  previous  steady-state  developments  of  the  Fourier- 
ray  method  are  extended  to  include  a  transient  source  of  mountain  waves.  The  method  is  illustrated  with 
an  initial  value  problem  in  which  the  background  flow  is  started  abruptly  from  rest  and  then  maintained  at 
steady  velocity.  The  resulting  wave  transience  is  modeled  in  a  simple  way.  All  rays  that  radiate  from  the 
mountain,  including  the  initial  rays,  are  assigned  the  full  amplitude  of  the  longtime  steady-state  solution. 
Time  dependence  comes  in  through  the  changing  position  of  the  initial  rays.  This  is  sufficient  to  account  for 
wave  transience  in  a  test  case,  as  demonstrated  by  comparison  with  simulations  from  a  mesoscale  numerical 
model. 


1.  Introduction 

Ray  solutions  are  often  expressed  in  spatial  coordi¬ 
nates.  Some  examples  for  mountain  waves  are  given  in 
Gjevik  and  Marthinsen  (1978),  Hines  (1988),  Shutts 
(1998),  and  Broad  (1999).  Ray  solutions  can  also  be 
expressed  in  Fourier-transform  coordinates  and  then 
mapped  into  a  spatial  solution  by  inverse  Fourier  trans¬ 
form.  For  waves  in  a  height-dependent  background,  the 
use  of  Fourier  transform  coordinates  simplifies  impor¬ 
tant  parts  of  the  ray  calculation,  including  the  ray  ini¬ 
tialization,  the  ray  tracing  itself,  and  the  correction  of 
caustics. 

We  believe  this  method  has  promise  for  mountain- 
wave  forecasting.  Three-dimensional  solutions  can  be 
calculated  in  minutes  on  a  standard  (1  GHz)  processor, 
and  at  much  higher  resolution  than  is  practical  with  a 
mesoscale  model.  The  method  handles  vertical  direc- 
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tional  wind  shear  in  the  presence  of  turning  points 
(Broutman  et  al.  2003)  and  critical  layers  (Broutman  et 
al.  2002).  However,  it  does  not  treat  finite- amplitude 
effects,  boundary  layer  effects,  flow  blocking,  or  hori¬ 
zontally  nonuniform  backgrounds. 

Previous  applications  of  the  method  have  been  re¬ 
stricted  to  longtime  steady-state  solutions  (Broutman  et 
al.  2002,  2003).  Here  we  develop  a  time-dependent  for¬ 
mulation  that  follows  the  evolution  of  a  transient  wave 
field.  This  can  lead  to  improved  forecasts  and  is  useful 
for  interpreting  observations  and  numerical  models. 
There  are  computational  benefits  as  well,  since  the 
transient  solution  can  sometimes  be  represented  on  a 
smaller  grid  and  does  not  have  the  resonant  singulari¬ 
ties  associated  with  the  trapped  waves  of  the  steady- 
state  solution  (see  section  2d). 

We  consider  an  initial  value  problem  in  which  the 
background  flow  is  started  abruptly  from  rest  and  then 
maintained  at  a  steady  velocity.  The  resulting  wave 
transience  is  modeled  in  a  simple  way.  All  rays  that 
emerge  from  the  mountain,  including  the  initial  rays, 
are  assigned  the  full  wave  amplitude  of  the  longtime 
steady-state  ray  solution.  The  time  dependence  comes 
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in  through  the  changing  position  of  the  initial  rays.  This 
approach  gives  a  reasonable  representation  of  the 
evolving  wave  field,  as  demonstated  in  section  3  by 
comparison  with  a  mesoscale  numerical  model  in  an 
idealized  case.  A  demonstration  for  a  more  realistic 
case  is  given  in  Eckermann  et  al.  (2006)  for  mountain 
waves  generated  by  the  island  of  Jan  Mayen. 

Our  theory  relates  to  a  transient  source  of  waves,  not 
a  transient  background  wind  or  stratification.  We  do 
envision  an  instant  start-up  of  the  wind  but  only  as  a 
means  to  generate  an  evolving  wave  field.  To  properly 
accommodate  a  transient  background  we  would  have  to 
allow  variations  in  the  wave  frequency  along  the  ray 
(e.g.,  Lott  and  Teitelbaum  1993),  and  we  would  have  to 
deal  with  a  more  general  range  of  caustic  conditions. 

We  refer  to  the  present  method  as  the  Fourier-ray 
method  to  indicate  that  the  spatial  solution  is  obtained 
by  a  Fourier  synthesis  of  ray  solutions.  This  is  the  third 
name  that  has  been  used  for  the  same  basic  approach  in 
three  papers.  The  reasons  for  the  name  changes  are 
discussed  at  the  end  of  section  2. 

2.  Theory 

We  modify  the  derivation  of  the  Fourier-ray  method 
in  Broutman  et  al.  (2003)  to  allow  for  wave  transience. 
We  consider  an  initial  value  problem  in  which  the  wind 
is  imposed  abruptly  at  time  t  =  0.  All  rays,  including  the 
initial  rays,  are  treated  as  though  their  wave  amplitudes 
instantly  satisfy  the  steady-state  solution  for  wave- 
action  conservation.  The  theory  does  allow  for  a 
gradual  ramping  up  of  the  wave  amplitude  to  a  steady- 
state  value,  rather  than  an  instantaneous  growth  to  the 
full  steady-state  value,  but  we  have  found  that  the  latter 
is  sufficient  to  obtain  a  good  comparison  with  numeri¬ 
cal  simulations. 

The  coordinate  system  is  x  =  (x,  y,  z)  with  z  vertical. 
The  background  wind  is  U  =  [U(z),  V(z),  0],  and  the 
buoyancy  frequency  is  N(z).  For  the  mountain  waves, 
k  =  ( k ,  /,  ±m)  is  the  wavenumber  vector,  and  w  = 
—kU  -  IV  is  the  intrinsic  frequency.  Our  notation  is 
that  m  >  0,  with  a  minus  sign  preceding  m  when  it  is 
appropriate  to  indicate  waves  with  downward  phase  ve¬ 
locity  and  upward  group  velocity.  The  gravity  wave  dis¬ 
persion  relation  is 

m2  =  (k2  +  l2)(N2/u>2  -  1).  (1) 

We  have  ignored  the  effects  of  compressibility  and  the 
earth’s  rotation  in  (1).  The  usual  density  scaling  of  the 
wave  amplitudes  can  be  incorporated  into  the  theory 
through  (8)  in  order  to  to  account  for  the  decrease  in 
atmospheric  density  with  height.  This  is  done  in  Eck¬ 


ermann  et  al.  (2006),  but  the  results  presented  in  sec¬ 
tion  3  below  are  calculated  in  the  Boussinesq  limit. 

We  express  the  ray  solution  in  Fourier-transform  co¬ 
ordinates,  with  the  Fourier  transform  taken  over  the 
horizontal  coordinates  only.  The  ray  solution  is  thus  a 
function  of  k ,  /,  z,  and  time  t.  Its  inverse  Fourier  trans¬ 
form  yields  an  appoximation  for  the  spatial  wave  solu¬ 
tion.  An  advantage  of  using  Fourier-transform  coordi¬ 
nates  for  the  ray  tracing  instead  of  the  more  usual  spa¬ 
tial  coordinates  is  that  k  and  l  are  constant  along  each 
ray  in  a  horizontally  uniform  background,  as  assumed 
here.  This  reduces  most  of  the  problem  to  a  one¬ 
dimensional  calculation  as  a  function  of  z.  There  is  also 
time  dependence,  which  comes  into  the  calculation  only 
in  tracking  the  height  of  the  initial  rays. 

We  work  here  in  terms  of  the  vertical  displacement, 
denoted  by  t)(x,  y,  z,  t)  in  spatial  coordinates  and  by 
r\(k,  /,  z,  t )  in  Fourier-transform  coordinates.  The 
former  is  the  inverse  Fourier  transform  of  the  latter: 

v(x,  y,z,t)=  1 1  v(k,  l,  z,  t)elikx+ly)  dk  dl.  (2) 

In  the  following  derivation,  we  treat  two  types  of  waves 
separately:  propagating  waves  without  turning  points, 
denoted  by  f)pr,  and  (vertically)  trapped  waves  with 
turning  points,  denoted  by  fjtr.  Thus, 

V  =  Vpr  +  V  tr  ■  (3) 

It  is  also  convenient  to  divide  the  trapped  waves  into 
contributions  from  upgoing  and  downgoing  rays: 

Vtr  Tlup  "b  ^Idn*  (^) 

a.  Propagating  waves 

The  solution  for  the  vertically  propagating  waves  has 
the  form 

rjpl(k,  l,  Z,  t)  =  FMG0/G)V2e‘4'.  (5) 

Except  for  the  factor  F,  which  is  an  amplitude  modu¬ 
lation  function  used  to  account  for  wave  transience,  this 
is  the  same  form  for  the  ray  solution  as  the  one  derived 
for  propagating  waves  in  Broutman  et  al.  (2003).  We 
now  define  the  terms  in  (5). 

The  wave  phase  in  (5)  is  given  by 

</>(£,  /,  z)  =  ~  I  m(k,  /,  z')  dz\  (6) 

J  o 

and  G  is  defined  by 

G  =  PN2cg3/w,  (7) 

with  G0  =  G(z  =  0).  The  vertical  group  velocity  is 
cg 3  =  da )ldm,  and  the  background  density  is  p(z). 
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The  term  (G0/G)1/2  in  (5)  ensures  that,  apart  from  the 
transient  effects  modeled  by  F,  the  waves  satisfy  the 
constancy  of  the  vertical  flux  of  wave  action 

pcg3A  =  constant,  (8) 

where  A  =  |  f)  1 2iV2/co  is  the  wave-action  density  per  unit 
mass.  This  constancy  of  the  vertical  flux  of  wave  action 
is  satisfied  only  in  the  Fourier  transform  domain.  The 
waves  in  the  spatial  domain  disperse  horizontally,  as 
well  as  vertically,  and  therefore  do  not  preserve  their 
vertical  flux  of  wave  action.  For  a  form  of  the  wave- 
action  solution  that  includes  the  effects  of  horizontal 
and  vertical  dispersion  (see  Shutts  1998;  Broad  1999). 

We  have  assumed  that  the  waves  satisfy  the  linear¬ 
ized  lower  boundary  condition 

fj(k,  l,z  =  0)  =  h(k,  /),  (9) 

where  h(k,  l )  is  the  Fourier  transform  of  the  topogra- 
phy  h(x,  y ). 

We  incorporate  wave  transience  through  the  function 
F(k ,  /,  z,  t )  in  (5).  This  is  the  only  quantity  on  the 
right-hand  side  of  (5)  that  depends  on  time.  It  can  be 
used  to  turn  the  wave  amplitude  on  or  off,  gradually  or 
instantly.  For  the  present  application,  we  turn  the  wave 
amplitude  on  instantly  and  permanently.  That  is,  we 
assume  that  the  waves  are  generated  from  an  initial 
state  of  rest  at  t  =  0,  and  we  assign  the  full  wave  am¬ 
plitude  of  the  steady-state  solution  to  every  ray.  Thus, 
we  set 


F(k,  /,  z,t)  =  H 


cg3(k ,  /,  z)  dt'  -  z 


(10) 


The  Heaviside  function  i/(£)  is  zero  for  £  <  0  and  unity 
for  £  >  0.  The  integral  in  the  argument  of  H  is  the 
height  at  time  t  of  the  initial  ray  (for  each  k ,  /)  gener¬ 
ated  at  the  mountain  at  t  =  0. 


b.  Trapped  waves 

Trapped  rays  continually  reflect  between  the  ground 
and  the  turning-point  height  zt,  where  6)  =  N.  For  each 
reflection,  a  new  term  must  be  added  to  the  ray  solu¬ 
tion.  Each  new  term  is  a  copy  of  the  previous  term,  but 
phase  shifted  by  the  amount 

0  =  2 4>-  tt/2.  (11) 


Here 


4)  =  I  mdz, 


is  the  change  in  phase  for  propagation  from  the  ground 
to  the  turning-point  height  zt(k,  /).  The  factor  of  —  7t/2 
in  (11)  is  the  sum  of  two  phase  shifts:  the  phase  shift  of 
+  7t/2  due  to  reflection  from  the  turning  point,  and  the 
phase  shift  of  —  77  due  to  reflection  from  the  ground. 
[The  turning-point  phase  shift  is  given  on  p.  397  of 
Lighthill  (1978)  and  in  section  10.1.3  of  Kravtsov  and 
Orlov(1999).] 

Thus,  after  the  initial  ray  crosses  z  on  its  way  up  M1 
times,  the  ray  solution  for  the  upgoing  rays  is 

Vup  =  fnG(/G)h-e'tSMi.  (13) 

The  sum  SM(k,  /,  z)  is  defined  by 


M-l 

+  tat) 

(14) 

n  =  0 

■y  _  iol) 

y  _  ’ 

(15) 

with  M  =  M1  for  the  upgoing  rays  and  M  =  M2  for  the 
downgoing  rays  (see  below).  When  M1  =  0,  that  is, 
before  the  initial  ray  crosses  z  for  the  first  time,  SMi  is 
defined  to  be  zero. 

The  positive  constant  a  in  the  expression  for  S  rep¬ 
resents  damping.  The  results  we  present  in  the  next 
section  are  in  viscid  (a  =  0),  but  we  have  included  this 
damping  term  to  use  in  an  intermediate  step  of  a  deri¬ 
vation  below  (section  2d),  and  because  it  might  be  use¬ 
ful  in  other  applications.  For  example,  Smith  et  al. 
(2002)  applied  the  same  form  of  damping  to  parameter¬ 
ize  the  partial  absorption  of  trapped  waves  by  a  stag¬ 
nant  lower  boundary  layer.  The  problem  was  reexam¬ 
ined  by  Smith  et  al.  (2006),  who  found  that  a  better 
description  of  the  process  should  include,  in  addition  to 
damping,  a  phase  shift  different  from  77  in  the  wave 
reflected  from  the  boundary  layer.  The  phase  shift  can 
be  incorporated  in  the  present  model  with  an  imaginary 
part  for  a. 

After  the  initial  ray  crosses  the  height  z  on  its  way 
down  M2  times,  the  ray  solution  for  the  downgoing  rays  is 

rjdn  =  ih{G0/G)1,2e^SM2.  (16) 

Here  cj)2  is  the  phase  of  the  initial  ray  after  its  first 
reflection  from  the  turning  point  but  before  its  first 
reflection  from  the  ground.  We  include  the  it  12  phase 
shift  from  the  first  turning-point  reflection  as  the  factor 
of  tat  the  start  of  the  right-hand  side  of  (16).  This  leaves 

4>2  =  <t>  ~  J  m(k,l,z')  dz'.  (17) 

The  sum  SM2  in  (16)  is  given  by  (14)  with  M  =  M2. 


0 


(12) 


2852 


MONTHLY  WEATHER  REVIEW 


Volume  134 


Wave  transience  thus  comes  in  through  the  crossing 
numbers  Mt  and  M2,  which  are  functions  of  k,  l  and 
which  increase  in  a  step-function  manner  with  time. 
Either  =  M2  +  1,  if  the  initial  ray  has  just  crossed  z 
on  its  way  up,  or  =  M2  if  the  initial  ray  has  just 
crossed  z  on  its  way  down. 

Our  most  general  form  for  the  ray  solution  for  the 
trapped  waves  is  then 


from  the  uniform  solution.  The  trapped  wave  solution 
can  then  be  written  as 

Vtr  =  VU  +  (M,  -  M2)h(G0/G)1/2e^+M^\  (23) 

where  r\u  is  given  by  (22).  When  Mx  =  M2  each  ray 
incident  on  the  turning  point  is  paired  with  a  ray  re¬ 
flected  from  the  turning  point.  Otherwise  the  second 
term  above  accounts  for  the  initial  ray. 


Vtr  =  h(Go/G)1/2(el+SMl  +  te'+*SM2). 


(18) 


d.  Steady-state  solutions 


c.  Caustic  correction  for  the  trapped  waves 

At  a  turning  point,  the  ray  solution  (18)  diverges 
because  cg3 ,  and  hence  G,  vanish.  This  is  a  caustic  sin¬ 
gularity.  The  simplest  caustics  are  correctable  with  an 
Airy  function.  One  can  either  match  the  ray  solution  to 
the  Airy  function  at  a  position  near  the  caustic  (e.g., 
section  4.11  of  Ligh thill  1978),  or  combine  the  ray  so¬ 
lution  and  the  Airy  function  into  a  uniform  represen¬ 
tation  that  is  valid  close  to,  and  far  from,  the  caustic 
(e.g.,  section  10.1.3  of  Kravtsov  and  Orlov  1999).  The 
latter  method  is  simple  to  implement,  so  we  will  use  it 
here. 

We  first  consider  one  pair  of  rays,  incident  upon  and 
reflected  from  the  caustic,  with  solutions  of  the  form  A 
exp(L(/>1)  and  iA  exp(u/>2),  respectively.  The  factor  of  l 
in  the  reflected  wave  solution  accounts  for  the  phase 
shift  of  7t/2  at  the  caustic.  The  uniform  solution  is  (sec¬ 
tion  10.1.3  of  Kravtsov  and  Orlov  1999) 

u  =  2m1/2(-r)1/4AAi(r)ei^~7T/4),  (19) 


where 


€  =  |  (4>i  +  4>2),  (20) 


r  = 


(21) 


The  Airy  function  argument  r  is  zero  at  the  caustic, 
negative  in  the  propagating  region  below  the  caustic, 
and  positive  in  the  evanescent  region  above  the  caustic. 

When  there  are  M2  downgoing  rays,  there  are  also 
M2  pairs  of  incident  and  reflected  rays.  We  can  accom¬ 
modate  M2  pairs  by  again  introducing  the  sum  SM2,  so 
that  the  uniform  solution  is 


f,u  =  2m1/2(-r)1/4AAi(r)e^/4)SM2,  (22) 

with  A  =  h(G0/G )1/2. 

If  M1  =  M2  +  1,  the  initial  ray  is  on  its  way  up  to  the 
turning  point  and  is  not  paired  with  a  reflected  ray  from 
the  turning  point.  It  thus  needs  to  be  counted  separately 


Our  time-dependent  solution  for  the  trapped  waves 
approaches  the  steady-state  solution  of  Broutman  et  al. 
(2003)  in  the  infinite-time  limit.  This  can  be  shown  by 
letting  M1  and  M2  approach  infinity  and  introducing  a 
small  amount  of  damping  in  (15),  through  the  positive 
constant  a,  so  that  the  sums  SMi  and  SMi  converge. 
Substituting  (15)  into  (18)  and  letting  a  — >  0  gives,  after 
some  algebra, 


.  /f  ,  ,  r~/ ^o\1/2  sin(4>  —  <fi  ~  tt/4)  _ 

7 ]tr(k,l,z,t  -^°°)  =  h[—  - 7 - .  (24) 

G  /  sin(c/>  -  7t/4) 


This  agrees  with  (18)  of  Broutman  et  al.  (2003).  For  the 
uniform  solution  (22),  the  limit  of  M1  and  M2  approach¬ 
ing  infinity  gives,  after  some  algebra, 


~(  G0 

huikfzd^00)  =  M  -q 


1/4  Ai(r) 
Ai(r0)  ’ 


(25) 


This  agrees  with  (A9)  of  Broutman  et  al.  (2003). 

Some  of  our  expressions  for  the  transient  solutions 
can  be  expressed  in  terms  of  a  time-dependent  factor 
times  the  steady-state  solution.  But  the  computational 
disadvantage  of  doing  this  is  that  the  steady-state  solu¬ 
tion  contains  resonant  singularities.  These  occur  for 
those  k ,  l  values  for  which  the  reflected  rays  are  per¬ 
fectly  in  phase  with  each  other.  They  correspond  to  the 
zeros  of  sin(4>  -  7t/4)  in  (24)  and  to  the  zeros  Ai(r0 )  in 
(25).  A  small  imaginary  component  added  to  the  wave- 
number  (or  frequency)  eliminates  this  singularity,  as 
was  done  in  Broutman  et  al.  (2003),  but  this  also  damps 
the  solution.  No  such  singularities  exist  in  the  time- 
dependent  solutions,  and  there  is  no  associated  damp¬ 
ing  in  our  solutions. 


e.  Terminology 

As  noted  in  the  introduction,  three  different  names 
have  been  applied  to  the  same  basic  method  used  here 
of  synthesizing  ray  solutions  by  inverse  Fourier  trans¬ 
form. 

In  Broutman  et  al.  (2002),  the  method  was  called 
Maslov’s  method,  after  a  Russian  mathematician  whose 
work  is  summarized  by  Maslov  and  Fedoriuk  (1981) 
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and  at  a  more  understandable  level  by  Brown  (2000). 
The  idea  is  that  as  the  rays  are  mapped  between  the 
spatial  and  Fourier-transform  domains,  the  locations  of 
caustics  change.  In  some  cases  the  caustics  in  one  do¬ 
main  disappear  entirely  in  the  other  domain.  This  oc¬ 
curs  in  the  hydrostatic  mountain-wave  model  studied 
by  Broutman  et  al.  (2002).  They  derived  the  ray  solu¬ 
tion  in  the  spatial  domain,  then  mapped  it  into  a  ray 
solution  in  the  Fourier  transform  domain  using  the  ray 
relation  between  position  and  wavenumber.  They  then 
mapped  back  to  the  spatial  domain  by  inverse  Fourier 
transform.  This  sequential  use  of  different  types  of 
mappings — first  by  the  ray  relations,  then  by  inverse 
Fourier  transform — is  why  one  finishes  with  a  different 
spatial  solution  than  the  original  spatial  ray  solution. 

This  is  the  approach  developed  originally  by  Maslov. 
But  the  application  of  Broutman  et  al.  (2002)  was 
simple:  there  were  no  caustics  in  the  Fourier-transform 
domain,  and  furthermore  the  problem  was  separable, 
and  thereby  reducible  to  a  one-dimensional  calculation 
in  the  Fourier-transform  domain.  Maslov  is  recognized 
for  his  work  on  the  more  complicated  case  of  nonsepa- 
rable  problems  with  caustics  in  both  the  spatial  and 
Fourier-transform  domains,  as  in  the  application  to  sur¬ 
face  gravity  waves  by  Brown  (2000). 

So  in  the  study  of  Broutman  et  al.  (2003),  the  name 
was  changed  to  a  “simplified  Fourier  method”  (see  also 
Broutman  and  Rottman  2004).  The  Fourier  method  in 
this  context  is  a  Fourier  synthesis  of  the  vertical  eigen¬ 
functions,  and  the  simplification  refers  to  that  fact  that 
the  vertical  eigenfunctions  were  approximated  by  ray 
theory.  In  the  present  paper,  we  changed  the  name  yet 
again  to  the  “Fourier-ray  method”  in  order  to  empha¬ 
size  the  involvement  of  ray  theory  and  its  enhanced  role 
here  in  predicting  the  evolving  positions  of  the  leading 
rays  for  the  time-dependent  solutions. 

3.  Results 

For  a  test  case,  we  consider  the  model  of  Wurtele  et 
al.  (1987,  hereafter  WSK).  The  stratification  is  uniform, 
and  the  mountain  waves  are  trapped  by  a  horizontal 
wind  that  increases  linearly  with  height.  We  have  al¬ 
ready  used  this  model  for  a  test  case  in  the  steady-state 
theory  of  Broutman  et  al.  (2003).  As  in  that  paper,  we 
use  a  nonlinear  mesoscale  model  for  verification.  It  has 
a  standard  second-order  finite-difference  discretization 
of  the  nonlinear  nonhydrostatic  equations  of  motion,  in 
either  the  anelastic  approximation  or,  as  used  here,  the 
Boussinesq  approximation.  The  model  is  essentially  the 
same  as  that  described  in  Lipps  and  Hemler  (1982). 

The  topography  is 

h(x,y)  =  h0a3/(x2  +  y2  +  a2)3/2,  (26) 


with  Fourier  transform 

h(k,l)  =  (h0a2/2TT)  exp [~a(k2  +  l2)1/2\  (27) 

The  mean  wind  is  given  by 

U(z)  =  U0(l  +  z/L ),  (28) 

where  U0  and  L  are  constants.  The  y-component  V  of 
the  wind  is  zero.  We  consider  the  tropospheric  model  in 
WSK’s  section  2a  with  the  parameter  values  of  WSK’s 
case  II.  These  are  N  =  0.01  s-1,  a  =  2.5  km,  h0  =  100  m, 
L  =  2.5  km,  and  U0  =  10  m  s-1. 

To  determine  and  M2,  the  time-dependent  num¬ 
bers  of  upgoing  and  downgoing  rays,  we  have  found  it 
adequate  to  use  the  following  approximation.  Each 
time  the  initial  ray  (generated  at  the  mountain  at  t  =  0) 
makes  one  complete  round  trip  from  the  ground  to  the 
turning  point  and  back  again,  we  increase  both  Mx  and 
M2  by  one.  We  use  a  rounded  estimate  given  by 

M1  =  M2  =  round(£/2T),  (29) 

where  T  is  the  propagation  time  from  the  ground  to  the 
turning  point,  and  t/2T  is  rounded  to  the  nearest  inte¬ 
ger. 

The  expression  for  T  can  be  derived  from  the  ray 
equation  dm/dt  =  ~kUz ,  where  d/dt  is  the  rate  of 
change  following  the  ray  (Lighthill  1978).  Using  (28), 
the  ray  equation  integrates  to  m  =  m0  -  kU0t/L.  Let¬ 
ting  m0  =  m(z  =  0),  and  identifying  m  =  0  as  the 
vertical  wavenumber  at  the  turning  point,  we  have 

T  =  m0L/kU0  (30) 

for  the  propagation  time  from  the  ground  to  the  turning 
point. 

We  show  solutions  for  the  vertical  velocity  w  = 
Dr\/Dt ,  where  D/Dt  =  d/dt  +  Ud/dx  is  the  advective  time 
derivative,  and  r\  is  the  vertical  displacement.  Since  w 
and  rj  are  related  by  w  =  —  icorj,  all  of  the  solutions  that 
we  have  derived  for  f)  need  only  to  be  multiplied  by 
—  ico  to  obtain  the  solution  h>.  An  inverse  Fourier  trans¬ 
form  of  w  gives  the  spatial  solution  w. 

For  the  Fourier-ray  model,  the  solution  is  computed 
by  a  discrete  Fourier-series  approximation  to  the  in¬ 
verse  Fourier  transform,  for  a  spatial  grid  of  512  points 
in  both  x  and  y.  (The  vertical  cross  sections  in  Fig.  2 
were  computed  at  a  1-km  vertical  grid  spacing.) 

For  the  numerical  model,  there  are  256  grid  points  in 
each  horizontal  direction,  and  32  grid  points  in  the  ver¬ 
tical  direction,  with  a  grid  spacing  of  1  km  in  all  direc¬ 
tions.  There  are  absorbing  sponges  near  the  top  and 
lateral  boundaries  of  the  numerical  model.  The  sponge 
near  the  top  boundary  is  based  on  a  Rayleigh  friction, 
with  the  form  Rqq'  in  the  governing  equations.  Here  q' 
is  a  perturbation  quantity  (velocity  or  potential  tern- 
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Fig.  1.  Comparison  of  (left)  the  Fourier-ray  solution  and  (right)  the  numerical  model  solution  for  the  vertical  velocity  w  in  the  WSK 
model,  after  1,  2,  and  3  h.  The  solutions  are  plotted  at  z  =  2.5  km.  Shaded  contours  are  positive,  the  zero  contour  is  omitted,  and  the 
contour  interval  is  0.025  ms-1. 


perature),  and  the  damping  coefficient  Rv  increases  lin¬ 
early  with  z  in  the  top  quarter  of  the  model  domain  to 
a  peak  value  of  0.004  s_1.  The  sponge  near  the  lateral 
boundaries  has  the  form  RhV2q' ,  where  Rh  increases 
linearly  across  the  final  10  horizontal  grid  points,  to  a 
peak  value  of  4500  m2  s-1. 

Figure  1  shows  a  comparison  of  (left  panels)  the  Fou¬ 
rier-ray  solution  and  (right  panels)  the  numerical  model 
solution  for  the  vertical  velocity  w  at  z  =  2.5  km  and  at 
times  t  =  1,  2,  and  3  h.  The  range  of  values  is  about 
-0.21  to  0.25  for  the  numerical  model,  and  about  -0.19 
to  0.23  for  the  Fourier-ray  method.  These  ranges  are 
established  by  t  =  1  h.  The  numerical  model  gives  peak 
values  that  are  a  few  percent  higher  than  those  of  the 
Fourier-ray  method,  even  though  the  Fourier-ray 
method  is  inviscid  and  the  numerical  model  has  some 
numerical  dissipation  (though  no  explicit  dissipation 
outside  the  sponge  regions).  We  are  not  sure  why  this  is 
the  case,  but  we  did  check  to  see  whether  evanescent 
modes  (with  co>  N  at  the  ground)  could  account  for  the 
difference  in  peak  values.  We  calculated  the  uniform 
solution  for  evanescent  modes  in  the  steady-state  limit, 
but  their  magnitude  was  insignificant  except  at  closer 
distances  to  the  mountain.  The  results  plotted  here  ex¬ 
clude  these  evanescent  modes. 

The  steady-state  Fourier-ray  solution  corresponding 


to  Fig.  1  here  was  given  in  Fig.  6  of  Broutman  et  al. 
(2003).  The  steady-state  solution  has  an  amplitude 
range  for  w  of  —0.18  to  0.21,  which  is  slightly  smaller 
than  the  range  predicted  here  by  the  transient  Fourier- 
ray  solution.  This  difference  is  a  result  of  the  damping 
used  in  the  steady-state  solution  to  remove  the  resonant 
singularity  for  the  trapped  waves. 

Figure  1  also  shows  good  agreement  of  the  predic¬ 
tions  for  the  horizontal  extent  of  the  wave  field.  For 
example,  at  t  =  2  h,  the  mountain  waves  in  both  the 
Fourier-ray  solution  and  the  numerical  solution  extend 
about  four  wavelengths  downwind  of  the  mountain  be¬ 
fore  starting  to  decay  rapidly  with  distance,  beyond 
about  x  =  60  km.  Near  this  leading  edge  of  the  wave 
field,  the  initial  rays  make  the  dominant  contribution  to 
the  solution,  and  the  Fourier-ray  method  is  likely  to  be 
inaccurate.  This  is  due  to  our  simple  procedure  for  giv¬ 
ing  the  initial  rays  the  full  wave  amplitude  of  the 
steady-state  solution,  and  also  to  our  use  of  the  round¬ 
ing  in  (29). 

At  t  =  1  h,  the  Fourier-ray  method  predicts  an- 
amolous  wave  amplitudes  in  the  region  upwind  of  the 
mountain  (x  <  0),  where  there  is  an  extra  leading  crest 
(shaded  contour)  that  does  not  occur  in  the  numerical 
solution.  The  reason  for  this  may  be  that  during  these 
early  stages  of  the  wave  field  development,  a  larger 


October  2006 


BROUTMAN  ET  AL. 


2855 


20 


x(km) 


x(km) 


120  140 


120  140 


fraction  of  the  rays  have  just  been  switched  on  to  full 
steady-state  amplitude  in  the  ray  tracing,  which  would 
cause  some  “ringing”  in  the  spatial  solution  near  the 
edges  of  the  wave  field.  As  time  increases,  more  rays 
are  correctly  represented  by  the  steady-state  ampli¬ 
tudes,  and  the  anamolous  upwind  crest  disappears,  as 
the  plots  indicate. 

Figure  2  shows  the  same  solutions  as  in  Fig.  1,  but  in 
vertical  cross  section  at  y  =  0.  Again  there  is  good 
agreement  between  the  Fourier-ray  solution  and  the 
numerical  model  solution. 

4.  Summary 

We  used  ray  tracing  in  a  Fourier-transform  domain 
to  calculate  the  wave  amplitudes  and  wave  phase,  and 
to  estimate  wave  propagation  times  in  order  to  account 
for  wave  transience.  These  ray  solutions  were  then  Fou¬ 
rier  synthesized  to  compute  the  spatial  solution.  We 
considered  an  initial  value  problem  in  which  the  flow  is 
started  abruptly  from  rest,  and  the  initial  rays  (and  all 
subsequent  rays)  are  assigned  the  full  wave  amplitude 
of  the  longtime  steady-state  solution.  This  gave  good 
agreement  with  the  predictions  of  a  nonlinear  meso- 
scale  model  in  an  example  taken  from  WSK.  A  case 
with  realisitic  topography,  winds  and  stratification  is 


considered  by  Eckermann  et  al.  (2006),  who  provide 
comparisons  of  the  transient  Fourier-ray  method  with 
numerical  simulations  and  satellite  observations. 
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